# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <stdlib.h>

#ifdef _OMP_PRESS_BC_
# include <omp.h>
#endif

extern double inlet_velocity(double r, double t, short index);
extern short interpol_1dim(const double *x, const double *u, unsigned short n, double xp, double *up);

# ifdef DENS_TDEP
extern double dens(double T);
# endif
# ifdef KVISC_TDEP
extern double kvisc(double T);
# endif
/*
# include <time.h>
# include <sys/timeb.h> */

extern void press_bc(const struct grid *g, short tflag, unsigned short iRK)
  {/* tflag = 0: current time step i.e. t; tflag = 1: next time-step i.e. t+dt. */
   short i, j, k, _k, nr, nz, nth, tid, nTRD;
   double ***u, ***v, ***w, ***T, *r, *th, *z, _dwdt_iRK, L, density0, t, kvisc0, _z[PRESS_INTRP_N], (*_p)[PRESS_INTRP_N];
// struct timeb t1, t2;

// ftime(&t1);
#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'press_bc()'...\n");
#  endif

   if( (iRK > 0 && tflag != 0) || (iRK == 0 && tflag != 1) || (iRK > 2) )
       {
        fprintf(g->lptr, "\nERROR: In 'press_bc()': (tflag, iRK) = (%d, %d).\nThis function is designed to handle only following cases:\n\t[tflag = 0, iRK = 1, 2] and [tflag = 1, iRK = 0]\n", tflag, iRK);
        exit(-1);
       }

   t = g->t + tflag*g->dt;

   nTRD = 1;
#  ifdef _OMP_PRESS_BC_
   nTRD = g->nTRD;
#  endif

/*** "T" is used in "discr.h-->kvisc() etc.", while discretizing momentum-equation terms. ***/
   u = g->u;   v = g->v;   w = g->w; T = g->T; r = g->r;   th = g->th;   z = g->z;  nr = g->nr;   nz = g->nz; nth = g->nth;
   density0 = g->density0; kvisc0 = g->kvisc0;

/* These are used in the interpolation of pressure at outlet. */
   _p = (double (*)[PRESS_INTRP_N])malloc(nTRD*PRESS_INTRP_N*sizeof(double));
   for(_z[0] = z[NzG], k = 1; k < PRESS_INTRP_N; k++) _z[k] = 0.5*(z[NzG-1+k]+z[NzG+k]);

#  ifdef _OMP_PRESS_BC_
   #pragma omp parallel private(i,j,k,_k,_dwdt_iRK,L,tid)
#  endif
     {
   /* At the INLET: Numan BC */
      k = NzG-1; /*** At the BOTTOM INTERFACE. It is used in evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < g->nr_nozzle + NrG; i++) for(j = NthG; j < nth + NthG; j++)
           {
            switch(iRK)
               {
                case 0:
                   _dwdt_iRK = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1);
                   break;

                case 1:
                   _dwdt_iRK = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                               + c2*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2);
                   break;

                case 2:
                   _dwdt_iRK = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                               + (a31+a32)*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2)\
                               + a32*c2*g->dt*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 3);
                   break;
               }

/*        if(iRK == 3) _dwdt_iRK = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                             + (a41+a42+a43)*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2)\
                             + (a42*C2+a43*a31+a43*a32)*g->dt*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 3)\
                             + a43*a32*C2*g->dt*g->dt*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 4); */

            L =  0 - _dwdt_iRK - ZMOMT_TERM_3 + ZMOMT_VISC - ZMOMT_TERM_5;
            g->p[i][j][NzG-1] = g->p[i][j][NzG] - density0*(z[NzG+1]-z[NzG])*L;
           }


   /* At the SOLID-BASE: Numan BC */
      k = NzG-1; /*** At the BOTTOM INTERFACE. It is used in evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG+g->nr_nozzle; i < NrG+nr-g->nr_gap; i++) for(j = NthG; j < nth+NthG; j++)
            {
             L = ZMOMT_VISC - ZMOMT_TERM_5;
             g->p[i][j][NzG-1] = g->p[i][j][NzG] - density0*(z[NzG+1]-z[NzG])*L;
            }

   /* At the OUTLET: Dirichlet B.C. : It can be used because mass flow rate at the outlet is not specified. Because
      the boundary condition specified at the outlet is Numan BC and Dirichlet B.C. at the outlet can be specified. */
#     ifdef _OMP_PRESS_BC_
      tid = omp_get_thread_num();
      #pragma omp for schedule(guided) nowait
#     else
      tid = 0;
#     endif
      for(i = NrG+nr-g->nr_gap; i < NrG+nr; i++) for(j = NthG; j < nth+NthG; j++)
            {
             _p[tid][0] = P_ATM + density0*GRAVITY*h_WATER;
             for(_k = 1; _k < PRESS_INTRP_N; _k++) _p[tid][_k] = g->p[i][j][NzG-1+_k];
             interpol_1dim(_z, _p[tid], PRESS_INTRP_N, 0.5*(z[NzG-1]+z[NzG]), &g->p[i][j][NzG-1]);
            }

   /* At the SIDE-WALL: Numan BC. */
      i = NrG+nr-1; /*** Used in the evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
            {
             L = RMOMT_VISC;
             g->p[NrG+nr][j][k] = g->p[NrG+nr-1][j][k] + density0*(r[NrG+nr]-r[NrG+nr-1])*L;
            }

   /* At the TOP-SURFACE. */
      k = NzG+nz-1; /*** Used in the evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided)
#     endif
      for(i = NrG; i < nr+NrG; i++) for(j = NthG; j < nth+NthG; j++)
            {
             L = ZMOMT_VISC - ZMOMT_TERM_5;
             g->p[i][j][NzG+nz] = g->p[i][j][NzG+nz-1] + density0*(z[NzG+nz]-z[NzG+nz-1])*L;
            }

   /* For the NEGATIVE radius: applying symmetric B.C. across the center. */
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided)
#     endif
      for(k = 0; k < nz+2*NzG; k++) for(j = NthG; j < NthG+nth; j++) for(i = 0; i < NrG; i++)
          g->p[NrG-1-i][j][k] = g->p[NrG+i][j][k];

   /* For ghost-cells in theta-direction: Applying symmetric BC. */
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided)
#     endif
      for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
         {
          g->p[i][NthG-j-1][k] = g->p[i][NthG+j][k];
          g->p[i][nth+NthG+j][k] = g->p[i][nth+NthG-1-j][k];
         }
     }

   free(_p);

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'press_bc()' ends.\n");
#  endif

   return;
  }
